%% This is the program to compute the price-sensitivity measures

%% Find MRR and Lambda
FIRMS=TORQSTOCK(:,2);
% Convert ISSM trades data into Matlab data with the following format:
%1	TKR	Decimal Format
%2	PERMNO	Permno
%3	DATE	YYYYMMDD
%4	SPM     Seconds past midnight
%5	PRICE	Dollars
%6	SIZE	Transaction size
%7	COND	Transaction cond
%8	OEXCH	Originating exchange
%9	PECXH	Present exchange

%% TT is raw data from ISSM trades
TT=csvread('ISSMTRADESTORQALL.csv');
%Convert 1,2 digit tickers with 3 digit
ID1=[6567
6576
6582
6665
6671
6776
6780
6785
6873
6880
7169
7269
7270
7383
7582
7767
7779
7788
7873
7884
8072
84
87
89
];
ID2=[656732
657632
658232
666532
667132
677632
678032
678532
687332
688032
716932
726932
727032
738332
758232
776732
777932
778832
787332
788432
807232
843232
873232
893232
];
for i=1:24
    clc
    i
    TT(find(TT(:,1)==ID1(i)),1)=ID2(i);
end

TTSORT=sortrows(TT,[1 4 3 2 5 6 7]);

%find PERMNOS
PERMNOS=NaN*ones(length(TTSORT),1);
for i=1:length(PERMNOS)
    clc
    i
    PERMNOS(i,1)=TORQSTOCK(find(TTSORT(i,1)==TORQSTOCK(:,1)),2);
end
DT=TTSORT(:,4).*10000+TTSORT(:,3).*100+TTSORT(:,2);
SPM=TTSORT(:,5)*60*60+TTSORT(:,6)*60+TTSORT(:,7);
TT2=[TTSORT(:,1) PERMNOS DT SPM TTSORT(:,8:12)];
TORQTRADES=sortrows(TT2,[1 3 4]);

%Convert ISSM trades data into Matlab data with the following
%format:
%1	TKR     Decimal Format
%2	PERMNO	Permno
%3	DATE	YYYYMMDD
%4	SPM     Seconds past midnight
%5	OFR     Offer price
%6	BID     Bid price
%7	MODE	Quote Condition
%8	OEXCH	Original exchange
%9	BIDSIZE	Bid size
%10	OFRSIZE	Offer size
%11	PEXCH	Present exchange

clear TT TTSORT TT2
%% TT is raw data from ISSM quotes
TT=csvread('ISSMQUOTESTORQALL.csv');
%Convert 1,2 digit tickers with 3 digit
ID1=[6567
6576
6582
6665
6671
6776
6780
6785
6873
6880
7169
7269
7270
7383
7582
7767
7779
7788
7873
7884
8072
84
87
89
];
ID2=[656732
657632
658232
666532
667132
677632
678032
678532
687332
688032
716932
726932
727032
738332
758232
776732
777932
778832
787332
788432
807232
843232
873232
893232
];
for i=1:24
    clc
    i
    TT(find(TT(:,1)==ID1(i)),1)=ID2(i);
end

TTSORT=sortrows(TT,[1 4 3 2 5 6 7]);

%find PERMNOS
PERMNOS=NaN*ones(length(TTSORT),1);
for i=1:length(PERMNOS)
    clc
    i
    PERMNOS(i,1)=TORQSTOCK(find(TTSORT(i,1)==TORQSTOCK(:,1)),2);
end

DT=TTSORT(:,4).*10000+TTSORT(:,3).*100+TTSORT(:,2);
SPM=TTSORT(:,5)*60*60+TTSORT(:,6)*60+TTSORT(:,7);
TT2=[TTSORT(:,1) PERMNOS DT SPM TTSORT(:,8:14)];
TORQQUOTES=sortrows(TT2,[1 3 4]);

%Placeholders
BSGHLAMBDA=NaN*ones(length(FIRMS),1);
MRRRHOGMM=NaN*ones(length(FIRMS),1);
MRRLAMBDAGMM=NaN*ones(length(FIRMS),1);
MRRTHETAGMM=NaN*ones(length(FIRMS),1);
MRRPHIGMM=NaN*ones(length(FIRMS),1);

%for i=1:length(FIRMS); %MRR GMM has a problem at
%i=45,49,50,80,81,94,100,116,118,141
for i=1:144
        clc
        i
    clear QID TID FQ FT TD QD ID1 BSDTEMP QDZ P1 ID2 XTEMP QTEMP VTEMP BSDTEMP2 PTEMP XT XT1 PC QT QT1 BSDT BSDT1 BSQT X Y
    QID=find(...
        (TORQQUOTES(:,8)==78)&...%NYSE quotes
        ((TORQQUOTES(:,7)==0)|(TORQQUOTES(:,7)==65)|(TORQQUOTES(:,7)==66)|(TORQQUOTES(:,7)==72)|(TORQQUOTES(:,7)==79)|(TORQQUOTES(:,7)==82))&... %BBO eligible
        (TORQQUOTES(:,5)>TORQQUOTES(:,6))&...   % Ask>Bid
        (TORQQUOTES(:,6)>0)&...                 %Positive prices
        (mod(TORQQUOTES(:,6),(1/16))==0)&...    %Prices in 16ths
        (mod(TORQQUOTES(:,5),(1/16))==0)&...   %Prices in 16ths
        (TORQQUOTES(:,2)==FIRMS(i))...         %Firm i
         );
    TID=find(...
        ((TORQTRADES(:,5)<=200)&(TORQTRADES(:,5)>=1))&... %Limits on max/min prices
        (TORQTRADES(:,8)==78)&... % NYSE Trades only
        (TORQTRADES(:,2)==FIRMS(i))&...  %firm i
        (TORQTRADES(:,7)==0)... % Regular trade condition codes
        );
    if length(TID)>50;  % Restrict to at least 50 observations for estimation
        FQ=TORQQUOTES(QID,3:6);
        FT=TORQTRADES(TID,3:6);
        TD=NaN*ones(length(FT),1); %MRR XT
        QD=NaN*ones(length(FT),1); %Huang Stoll Q
        for t=1:length(TD);
            clear ID1;
            ID1=max(find((FQ(:,1)==FT(t,1))&(FQ(:,2)<FT(t,2)-5))); %Most recent quote on same day with 5s lapse
            if isempty(ID1)==0;
                if FT(t,3)==FQ(ID1,3); %If trade price = Ask, Trade direction =1
                    TD(t,1)=1;
                elseif FT(t,3)==FQ(ID1,4); %If trade price = bid, Trade direction =-1
                    TD(t,1)=-1;
                elseif (FT(t,3)<FQ(ID1,3)) && (FT(t,3)>FQ(ID1,4)); % If trade price is inside bid/ask, direction = 0
                    TD(t,1)=0;
                end
                if FT(t,3)>mean(FQ(ID1,3:4)); %If trade price > midpoint, Trade direction =1
                    QD(t,1)=1;
                elseif FT(t,3)<mean(FQ(ID1,3:4)); %If trade price < midpoint, Trade direction =-1
                    QD(t,1)=-1;
                elseif FT(t,3)==mean(FQ(ID1,3:4)); % If trade price is equal to midpoint, direction = 0
                    QD(t,1)=0;
                end
            end
        end
        BSDTEMP=QD; %Brennan Subrahmanyam q (includes tick test)
        clear QDZ IDQ
        QDZ=find(QD==0);
        for z=1:length(QDZ)
            clear IDQ
            IDQ=find(FT(1:QDZ(z)-1,3)~=FT(QDZ(z),3),1,'last');
            if FT(IDQ,3)<FT(QDZ(z),3)
                BSDTEMP(QDZ(z))=1;
            elseif FT(IDQ,3)>FT(QDZ(z),3);
                BSDTEMP(QDZ(z))=-1;
            end
        end
        clear P1 ID2 XTEMP PTEMP XT XT1 PC QTEMP VTEMP 
        P1=[NaN;FT(2:end,3)-FT(1:end-1,3)]; %Price changes (2:end of Firm trades)
        ID2=1+find((FT(2:end,1)-FT(1:end-1,1)==0)&(isnan(TD(2:end))==0)); % Not Overnight trades
        XTEMP=TD(ID2); % Order direction without overnight trades & no quotes outside spread
        BSDTEMP2=BSDTEMP(ID2); % Brennan Subra D
        PTEMP=P1(ID2); % Price changes without overnight trades & no quotes outside spread
        VTEMP=FT(ID2,4); %Volume without overnight trades & no quotes outside spread
        XT =XTEMP(2:end);   %Order direction at t MRR
        XT1=XTEMP(1:end-1); %Lagged order direction (t-1)
        PC =PTEMP(2:end);   %Price change at t
        BSDT=BSDTEMP2(2:end);
        BSDT1=BSDTEMP2(1:end-1);
        BSQT=BSDTEMP2(2:end).*VTEMP(2:end);
        %MRR 97 Estimates
        clear X Y MRROLS1 MRROLS2 MRROLS4 MRROLSB
        Y=[PC 1-abs(XT) XT.*XT1];
        X=[ones(length(XT),1) XT -XT1 XT.^2];
        if (i~=45)&(i~=49)&(i~=50)&(i~=80)&(i~=81)&(i~=94)&(i~=100)&(i~=116)&(i~=118)&(i~=141)         %i=45,49,50,80,81,94,100,116,118,141 have problems with GMM...?
            clear gmmopt MRR
            gmmopt.infoz.momt='mrrgmmm';  % residual and moment function
            gmmopt.prt='0'; %suppress printing
            gmmopt.plot='0'; %suppress plots
            MRR = gmm(zeros(5,1),gmmopt,Y,X,X);
            %b=[    alpha(constant) 
            %       phi + theta
            %       phi + rho * theta
            %       lambda
            %       rho]
            MRRRHOGMM(i,1)=MRR.b(5);
            MRRTHETAGMM(i,1)=(MRR.b(2)-MRR.b(3))/(1-MRRRHOGMM(i,1));
            MRRPHIGMM(i,1)=MRR.b(2)-MRRTHETAGMM(i,1);
        end
        %Brennan Subra 95 lambda estimates
        % Glosten Harris Spec.
        clear X Y
        Y=PC;
        X=[ones(length(PC),1) BSQT BSDT-BSDT1];
        clear gmmopt BSGH
        gmmopt.infoz.momt='lingmmm';  % residual and moment function
        gmmopt.prt='0'; %suppress printing
        gmmopt.plot='0'; %suppress plots
        BSGH = gmm(zeros(3,1),gmmopt,Y,X,X);
        %b=[S/2
        %   lambda*S/2]
        BSGHLAMBDA(i,1)=BSGH.b(2);
     end % 
end % End firm loop

MRRRGMM=MRRTHETAGMM./(MRRTHETAGMM+MRRPHIGMM); %Fraction of implied spread due to asymmetric info
TORQSTOCK(:,10:11)=[MRRRGMM BSGHLAMBDA];
